/*
 * Nbody.java
 *
 * Created on December 3, 2007, 4:29 AM
 *
 * To change this template, choose Tools | Template Manager
 * and open the template in the editor.
 */

package simfinal;

import com.jme.math.Vector3f;

/**
 *
 * @author jkeister
 */
public class Nbody 
{
    Particle[] system;
    double eps;
    double energy;
    double time, tnext;
    
    /** Creates a new instance of Nbody */
    public Nbody(Particle[] x, double e) 
    {
        
        System.out.println("Inside Nbody");
        
        system = x;
        eps = e;
        time = 0;
        tnext = 0;
        
        setup();
    }
    
    private void setup()
    {
        int N = system.length;
        // set up force and first derivative for each body
        for(int i=0; i<N; ++i)
        {
            Vector a, adot;
            for(int j=0; j<N; ++j)
            {
                if(i==j) continue;

                a = new Vector(system[j].x0);
                a.sub(system[i].x0);
                adot = new Vector(system[j].x0dot);
                adot.sub(system[i].x0dot);
                double r = 1.0/(a.dot(a)+eps);
                double a8 = system[j].mass * r *Math.sqrt(r);
                double a9 = 3.0*r*adot.dot(a);
                system[i].f.addm(a, a8);
                system[i].fdot.addm(adot, a8);
                system[i].fdot.addm(a, -a9*a8);
            }
        }

        // set up 2nd and 3rd force derivatives
        for(int i=0; i<N; ++i)
        {
            Vector a, adot, f, fdot;
            for(int j=0; j<N; ++j)
            {
                if(i==j) continue;

                a = new Vector(system[j].x0);
                a.sub(system[i].x0);
                adot = new Vector(system[j].x0dot);
                adot.sub(system[i].x0dot);
                f = new Vector(system[j].f);
                f.sub(system[i].f);
                fdot = new Vector(system[j].fdot);
                fdot.sub(system[i].fdot);

                double r = 1.0/(a.dot(a)+eps);
                double a14 = system[j].mass * r *Math.sqrt(r);
                double a15 = r*adot.dot(a);
                double a16 = (adot.dot(adot) + a.dot(f))*r + a15*a15;
                double a17 = (9.0*adot.dot(f) + 3.0*a.dot(fdot))*r +
                    a15*(9.0*a16 - 12*a15*a15);

                Vector f1dot = new Vector(adot);
                f1dot.addm(a, -3.0*a15);

                Vector f2dot = new Vector(f);
                f2dot.addm(f1dot, -6.0*a15);
                f2dot.addm(a, -3.0*a16);
                f2dot.mult(a14);
                system[i].d2.add(f2dot);

                Vector f3dot = new Vector(fdot);
                f3dot.addm(f1dot, -9.0*a16);
                f3dot.addm(a, -a17);
                f3dot.mult(a14);
                system[i].d3.add(f3dot);
                system[i].d3.addm(f2dot, -9.0*a15);
            } // j
            system[i].setStep();
        }
    }
    
    public Vector3f getPosition(int index)
    {
        return system[index].x0.toVector3f();
    }

    public Particle[] advance(double deltat)
   {
       // compute initial energy
       energy = 0.0;
       int N = system.length;
       int i, j;
       for(i=0; i<N; ++i)
       {
            energy += system[i].kinetic(tnext);
            for(j=0; j<i; ++j)
            {
                Vector s = new Vector(system[i].x);
                s.sub(system[j].x);
                double pe = system[i].mass*system[j].mass/Math.sqrt(s.dot(s)+eps);
                energy -= pe;
                system[i].energy -= pe/2.0;
                system[j].energy -= pe/2.0;
            }
       }
       tnext += deltat;
       boolean end = false;

       do{
            // find body with shortest timescale
            time = 1.0e38;
            for(j=0; j<N; ++j)
            {
                double t = system[j].timescale();
                if(time < t) continue;
                i = j;
                time = t;
            }
            
            if(time > tnext)
            {
                time = tnext;
                end = true;
            }

            // predict coordinates to 1st order in force
            for(j=0; j<N; ++j)
            {
                system[j].firstorder(time);
            }

            // set the critical body accurately
            system[i].thirdorder(time);

            // get the force on this body
            Vector f1 = new Vector();
            for(j=0; j<N; ++j)
            {
                if(j==i) continue;
                Vector a = new Vector(system[j].x);
                a.sub(system[i].x);
                double r = 1.0/(a.dot(a)+eps);
                double c = system[j].mass*r*Math.sqrt(r);
                f1.addm(a, c);
            }

            // and update all associated quantities
            system[i].update(time, f1);

        }while(!end);
        return system;
    }
    
}
